%% Estimation of parameters from noisy data

%1st order polynomial
a=1;
b=5;

x=(0:100);
y=a+b*x;

%Gaussian (normal) noise
sigma=20;
y_N=normrnd(y,sigma);

figure
plot(x,y)
xlabel('x')
ylabel('y')
hold on
plot(x,y_N,'r.')

legend('Known Model','Noisy Data')

%% Estimate from noisy data

%Use fminsearch to minimize negative log likelihood

%Here is a good guess for the intial search values. 
Theta0=[min(y_N); (y_N(end)-y_N(1))/(x(end)-x(1))] 

%Check that these values make a model that isn't too far off
% See 'firstOPoly_Norm'
[~,Model] = firstOPoly_Norm(Theta0,x,y_N)

figure
plot(x,y_N,'r.')
xlabel('x')
ylabel('y')
hold on
plot(x,Model)
legend('Noisy Data','intial Model')

%Setup fminsearch and find optimum values. 
Options=optimset();
[Theta_Found]=fminsearch(@firstOPoly_Norm,Theta0,Options,x,y_N)

%Get the model generated from found parameters
[~,Model_found] = firstOPoly_Norm(Theta_Found,x,y_N)

%Plot data and models.  
figure
plot(x,y_N,'r.')
xlabel('x')
ylabel('y')
hold on
plot(x,Model)
plot(x,Model_found,'g')
legend('Noisy Data','intial Model','found model')

%% Repeat for Poisson data

%1st order polynomial
a=1;
b=5;

x=(0:100);
y=a+b*x;

%Noisy data
y_N=poissrnd(y);

%Intial guess
Theta0=[min(y_N); (y_N(end)-y_N(1))/(x(end)-x(1))] 
% See 'firstOPoly_Poiss'
[~,Model] = firstOPoly_Poiss(Theta0,x,y_N)

%Setup fminsearch and find optimum values. 
Options=optimset();
[Theta_Found]=fminsearch(@firstOPoly_Poiss,Theta0,Options,x,y_N)

%Get the model generated from found parameters
[~,Model_found] = firstOPoly_Poiss(Theta_Found,x,y_N)

%Plot data and models.  
figure
plot(x,y_N,'r.')
xlabel('x')
ylabel('y')
hold on
plot(x,Model)
plot(x,Model_found,'g')
legend('Noisy Data','intial Model','found model')



